Licensed under the Creative Commons attribution-noncommercial license, http://creativecommons.org/licenses/by-nc/3.0/. Please share and remix noncommercially, mentioning its origin.
CC-BY_NC

Produced in R version 3.2.4 using pomp version 1.4.7.1.

Objectives

  1. To demonstrate the use of diagnostic probes for model criticism
  2. To teach some forecasting methods based on POMP models

These objectives will be achieved using a recent study (King et al. 2015), all codes for which are available on datadryad.org.

Model and data

An emerging infectious disease outbreak

Let’s situate ourselves at the beginning of October 2014. The WHO situation report contained data on the number of cases in each of Guinea, Sierra Leone, and Liberia. Key questions included:

  1. How fast will the outbreak unfold?
  2. How large will it ultimately prove?
  3. What interventions will be most effective?

Download the data from the WHO Situation Report of 1 October 2014:

base_url <- "http://kingaa.github.io/sbied/"
read.csv(paste0(base_url,"data/ebola_data.csv"),stringsAsFactors=FALSE,
         colClasses=c(date="Date")) -> dat
sapply(dat,class)
##        week        date     country       cases 
##   "integer"      "Date" "character"   "numeric"
head(dat)
##   week       date country cases
## 1    1 2014-01-05  Guinea 2.244
## 2    2 2014-01-12  Guinea 2.244
## 3    3 2014-01-19  Guinea 0.073
## 4    4 2014-01-26  Guinea 5.717
## 5    5 2014-02-02  Guinea 3.954
## 6    6 2014-02-09  Guinea 5.444

Supplementing these data are population estimates for the three countries.

## Population sizes in Guinea, Liberia, and Sierra Leone (census 2014)
populations <- c(Guinea=10628972,Liberia=4092310,SierraLeone=6190280)
dat %>%
  ggplot(aes(x=date,y=cases,group=country,color=country))+
  geom_line()

An SEIR model with gamma-distributed latent and infectious periods

Many of the early modeling efforts used variants on the simple SEIR model. Here, we’ll focus on a variant that attempts a more accurate description of the duration of the latent period. Specifically, this model assumes that the amount of time an infection remains latent is \[\mathrm{LP} \sim {\mathrm{Gamma}\left(m,\frac{1}{m\,\alpha}\right)},\] where \(m\) is an integer. This means that the latent period has expectation \(1/\alpha\) and variance \(1/(m\,\alpha)\). In this document, we’ll fix \(m=3\).

We implement Gamma distributions using the so-called linear chain trick.

Process model simulator

rSim <- Csnippet('
  double lambda, beta;
  double *E = &E1;
  beta = R0 * gamma; // Transmission rate
  lambda = beta * I / N; // Force of infection
  int i;

  // Transitions
  // From class S
  double transS = rbinom(S, 1.0 - exp(- lambda * dt)); // No of infections
  // From class E
  double transE[nstageE]; // No of transitions between classes E
  for(i = 0; i < nstageE; i++){
    transE[i] = rbinom(E[i], 1.0 - exp(-nstageE * alpha * dt));
  }
  // From class I
  double transI = rbinom(I, 1.0 - exp(-gamma * dt)); // No of transitions I->R

  // Balance the equations
  S -= transS;
  E[0] += transS - transE[0];
  for(i=1; i < nstageE; i++) {
    E[i] += transE[i-1] - transE[i];
  }
  I += transE[nstageE-1] - transI;
  R += transI;
  N_EI += transE[nstageE-1]; // No of transitions from E to I
  N_IR += transI; // No of transitions from I to R
')

Deterministic skeleton

The deterministic skeleton is an ODE.

skel <- Csnippet('
  double lambda, beta;
  const double *E = &E1;
  double *DE = &DE1;
  beta = R0 * gamma; // Transmission rate
  lambda = beta * I / N; // Force of infection
  int i;

  // Balance the equations
  DS = - lambda * S;
  DE[0] = lambda * S - nstageE * alpha * E[0];
  for (i=1; i < nstageE; i++)
    DE[i] = nstageE * alpha * (E[i-1]-E[i]);
  DI = nstageE * alpha * E[nstageE-1] - gamma * I;
  DR = gamma * I;
  DN_EI = nstageE * alpha * E[nstageE-1];
  DN_IR = gamma * I;
')

Measurement model: overdispersed count data

\(C_t | H_t\) is negative binomial with \({\mathbb{E}\left[{C_t|H_t}\right]} = \rho\,H_t\) and \({\mathrm{Var}\left[{C_t|H_t}\right]} = \rho\,H_t\,(1+k\,\rho\,H_t)\).

dObs <- Csnippet('
  double f;
  if (k > 0.0)
    f = dnbinom_mu(nearbyint(cases),1.0/k,rho*N_EI,1);
  else
    f = dpois(nearbyint(cases),rho*N_EI,1);
  lik = (give_log) ? f : exp(f);
')

rObs <- Csnippet('
  if (k > 0) {
    cases = rnbinom_mu(1.0/k,rho*N_EI);
  } else {
    cases = rpois(rho*N_EI);
  }')

Parameter transformations

toEst <- Csnippet('
  const double *IC = &S_0;
  double *TIC = &TS_0;
  TR0 = log(R0);
  Trho = logit(rho);
  Tk = log(k);
  to_log_barycentric(TIC,IC,4);
')

fromEst <- Csnippet('
  const double *IC = &S_0;
  double *TIC = &TS_0;
  TR0 = exp(R0);
  Trho = expit(rho);
  Tk = exp(k);
  from_log_barycentric(TIC,IC,4);
')

The following function constructs a pomp object to hold the data for any one of the countries.

ebolaModel <- function (country=c("Guinea", "SierraLeone", "Liberia"),
                        timestep = 0.1, nstageE = 3) {

  ctry <- match.arg(country)
  pop <- unname(populations[ctry])
  nstageE <- as.integer(nstageE)

  globs <- paste0("static int nstageE = ",nstageE,";")

  dat <- subset(dat,country==ctry,select=-country)

  ## Create the pomp object
  dat %>% 
    extract(c("week","cases")) %>%
    pomp(
      times="week",
      t0=min(dat$week)-1,
      globals=globs,
      statenames=c("S","E1","I","R","N_EI","N_IR"),
      zeronames=c("N_EI","N_IR"),
      paramnames=c("N","R0","alpha","gamma","rho","k",
                   "S_0","E_0","I_0","R_0"),
      nstageE=nstageE,
      dmeasure=dObs, rmeasure=rObs,
      rprocess=discrete.time.sim(step.fun=rSim, delta.t=timestep),
      skeleton=skel, skeleton.type="vectorfield",
      toEstimationScale=toEst,
      fromEstimationScale=fromEst,
      initializer=function (params, t0, nstageE, ...) {
        all.state.names <- c("S",paste0("E",1:nstageE),"I","R","N_EI","N_IR")
        comp.names <- c("S",paste0("E",1:nstageE),"I","R")
        x0 <- setNames(numeric(length(all.state.names)),all.state.names)
        frac <- c(params["S_0"],rep(params["E_0"]/nstageE,nstageE),params["I_0"],params["R_0"])
        x0[comp.names] <- round(params["N"]*frac/sum(frac))
        x0
      }
    ) -> po
}

ebolaModel("Guinea") -> gin
ebolaModel("SierraLeone") -> sle
ebolaModel("Liberia") -> lbr

Parameter estimates

King et al. (2015) estimated parameters for this model for each country. A large Latin hypercube design was used to initiate a large number of iterated filtering runs. Profile likelihoods were computed for each country against the parameters \(k\) (the measurement model overdispersion) and \(R_0\) (the basic reproductive ratio). Full details are given on the datadryad.org site. The following loads the results of these calculations.

options(stringsAsFactors=FALSE)
profs <- read.csv(paste0(base_url,"/ebola/ebola-profiles.csv"))

The following plots the profile likelihoods. The horizontal line represents the critical value of the likelihood ratio test for \(p=0.01\).

library(reshape2)
library(plyr)
library(magrittr)
library(ggplot2)
theme_set(theme_bw())

profs %>% 
  melt(id=c("profile","country","loglik")) %>%
  subset(variable==profile) %>%
  ddply(~country,mutate,dll=loglik-max(loglik)) %>%
  ddply(~country+profile+value,subset,loglik==max(loglik)) %>% 
  ggplot(mapping=aes(x=value,y=dll))+
  geom_point(color='red')+
  geom_hline(yintercept=-0.5*qchisq(p=0.99,df=1))+
  facet_grid(country~profile,scales='free')+
  labs(y=expression(l))

Diagnostics

Parameter estimation is the process of finding the parameters that are “best”, in some sense, for a given model, from among the set of those that make sense for that model. Model selection, likewise, aims at identifying the “best” model, in some sense, from among a set of candidates. One can do both of these things more or less well, but no matter how carefully they are done, the best of a bad set of models is still bad.

Lets’ investigate the model here, at its maximum-likelihood parameters, to see if we can identify problems. The guiding principle in this is that, if the model is “good”, then the data are a plausible realization of that model. Therefore, we can compare the data directly against model simulations. Moreover, we can quantify the agreement between simulations and data in any way we like. Any statistic, or set of statistics, that can be applied to the data can also be applied to simulations. Shortcomings of the model should manifest themselves as discrepancies between the model-predicted distribution of such statistics and their value on the data.

pomp provides tools to facilitate this process. Specifically, the probe function applies a set of user-specified probes or summary statistics, to the model and the data, and quantifies the degree of disagreement in several ways.

Let’s see how this is done using the model for the Guinean outbreak.

library(pomp)
library(plyr)
library(reshape2)
library(magrittr)
options(stringsAsFactors=FALSE)

profs %>%
  subset(country=="Guinea") %>%
  subset(loglik==max(loglik),
         select=-c(loglik,loglik.se,country,profile)) %>%
  unlist() -> coef(gin)

simulate(gin,nsim=20,as.data.frame=TRUE,include.data=TRUE) %>% 
  mutate(date=min(dat$date)+7*(time-1),
         is.data=ifelse(sim=="data","yes","no")) %>% 
  ggplot(aes(x=date,y=cases,group=sim,color=is.data,
         alpha=is.data))+
  geom_line()+
  guides(color=FALSE,alpha=FALSE)+
  scale_color_manual(values=c(no=gray(0.6),yes='red'))+
  scale_alpha_manual(values=c(no=0.5,yes=1))

The simulations appear to be growing a bit more quickly than the data. Let’s try to quantify this. First, we’ll write a function that estimates the exponential growth rate by linear regression. Then, we’ll apply it to the data and to 500 simulations.

growth.rate <- function (y) {
  cases <- y["cases",]
  fit <- lm(log1p(cases)~seq_along(cases))
  unname(coef(fit)[2])
}
probe(gin,probes=list(r=growth.rate),nsim=500) %>% plot()

Do these results bear out our suspicion that the model and data differ in terms of growth rate?

The simulations also appear to be more highly variable around the trend than do the data.

growth.rate.plus <- function (y) {
  cases <- y["cases",]
  fit <- lm(log1p(cases)~seq_along(cases))
  c(r=unname(coef(fit)[2]),sd=sd(residuals(fit)))
}
probe(gin,probes=list(growth.rate.plus),
      nsim=500) %>% plot()

Let’s also look more carefully at the distribution of values about the trend using the 1st and 3rd quantiles. Also, it looks like the data are less jagged than the simulations. We can quantify this using the autocorrelation function (ACF).

log1p.detrend <- function (y) {
  cases <- y["cases",]
  y["cases",] <- as.numeric(residuals(lm(log1p(cases)~seq_along(cases))))
  y
}

probe(gin,probes=list(
  growth.rate.plus,
  probe.quantile(var="cases",prob=c(0.25,0.75)),
  probe.acf(var="cases",lags=c(1,2,3),type="correlation",
            transform=log1p.detrend)
),nsim=500) %>% plot()

Exercise: the SEIR model for the Sierra Leone outbreak

Apply probes to investigate the extent to which the model is an adequate description of the data from the Sierra Leone outbreak. Have a look at the probes provided with pomp: ?basic.probes. Try also to come up with some informative probes of your own. Discuss the implications of your findings.

Forecasting

Up to now, we’ve primarily focused on using POMP models to answer scientific questions. Of course, we can also use them to make forecasts. The key issues are to do with quantifying the forecast uncertainty. This arises from four sources:

  1. measurement error
  2. process noise
  3. parametric uncertainty
  4. structural uncertainty

Here, we’ll explore how we can account for the first three of these in making forecasts for the Sierra Leone outbreak.

library(pomp)
library(plyr)
library(reshape2)
library(magrittr)
options(stringsAsFactors=FALSE)

set.seed(988077383L)

## forecast horizon
horizon <- 13

profs %>%
  subset(country=="SierraLeone") %>%
  subset(loglik==max(loglik),
         select=-c(loglik,loglik.se,country,profile)) %>%
  unlist() -> mle

## Weighted quantile function
wquant <- function (x, weights, probs = c(0.025,0.5,0.975)) {
  idx <- order(x)
  x <- x[idx]
  weights <- weights[idx]
  w <- cumsum(weights)/sum(weights)
  rval <- approx(w,x,probs,rule=1)
  rval$y
}

profs %>% 
  subset(country=="SierraLeone",
         select=-c(country,profile,loglik.se)) %>%
  subset(loglik>max(loglik)-0.5*qchisq(df=1,p=0.99)) %>%
  melt(variable.name="parameter") %>%
  ddply(~parameter,summarize,
        min=min(value),max=max(value)) %>%
  subset(parameter!="loglik") %>%
  melt(measure=c("min","max")) %>%
  acast(parameter~variable) -> ranges

params <- sobolDesign(lower=ranges[,'min'],
                      upper=ranges[,'max'],
                      nseq=20)
plot(params)

library(foreach)
library(doParallel)
library(iterators)

registerDoParallel()

set.seed(887851050L,kind="L'Ecuyer")

foreach(p=iter(params,by='row'),
        .inorder=FALSE,
        .combine=rbind,
        .options.multicore=list(preschedule=TRUE,set.seed=TRUE)
) %dopar%
{
  library(pomp)
  
  M1 <- ebolaModel("SierraLeone")
  pf <- pfilter(M1,params=unlist(p),Np=2000,save.states=TRUE)
  pf$saved.states %>% tail(1) %>% melt() %>% 
    dcast(rep~variable,value.var="value") %>%
    ddply(~rep,summarize,S_0=S,E_0=E1+E2+E3,I_0=I,R_0=R) %>%
    melt(id="rep") %>% acast(variable~rep) -> x
  
  pp <- parmat(unlist(p),ncol(x))
  
  simulate(M1,params=pp,obs=TRUE) %>%
    melt() %>%
    mutate(time=time(M1)[time],
           period="calibration",
           loglik=logLik(pf)) -> calib

    M2 <- M1
  time(M2) <- max(time(M1))+seq_len(horizon)
  timezero(M2) <- max(time(M1))
  
  pp[rownames(x),] <- x
  
  simulate(M2,params=pp,obs=TRUE) %>%
    melt() %>%
    mutate(time=time(M2)[time],
           period="projection",
           loglik=logLik(pf)) -> proj
  
  rbind(calib,proj)
} %>% subset(variable=="cases",select=-variable) %>%
  mutate(weight=exp(loglik-mean(loglik))) %>%
  arrange(time,rep) -> sims

ess <- with(subset(sims,time==max(time)),weight/sum(weight))
ess <- 1/sum(ess^2); ess
## [1] 10282.3
sims %>% ddply(~time+period,summarize,prob=c(0.025,0.5,0.975),
               quantile=wquant(value,weights=weight,probs=prob)) %>%
  mutate(prob=mapvalues(prob,from=c(0.025,0.5,0.975),
                        to=c("lower","median","upper"))) %>%
  dcast(period+time~prob,value.var='quantile') %>%
  mutate(date=min(dat$date)+7*(time-1)) -> simq
simq %>% ggplot(aes(x=date))+
  geom_ribbon(aes(ymin=lower,ymax=upper,fill=period),alpha=0.3,color=NA)+
  geom_line(aes(y=median,color=period))+
  geom_point(data=subset(dat,country=="SierraLeone"),
             mapping=aes(x=date,y=cases),color='black')+
  labs(y="cases")

Back to course homepage

R codes for this document

References

King, A. A., M. Domenech de Cellès, F. M. G. Magpantay, and P. Rohani. 2015. Avoidable errors in the modelling of outbreaks of emerging pathogens, with special reference to ebola. Proc. R. Soc. Lond. B 282:20150347.